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NCC 2-585 
ANNUAL REPORT 

Technical Proposal 

The objective of the present work is to develop, verify and incor- 
porate two-equation turbulence models which account for the effects of 
compressibility at high speeds into a three-dimensional Reynolds av- 
eraged Navier-Stokes (RANS) code and to provide documented model 
descriptions and numerical procedures so that they can be implemented 
into the NASP CFD codes. 

Computational Fluid Dynamics (CFD) is recognized as a signif- 
icant engineering design tool in modern hypersonic projects such as 
the National Aerospace Plane (NASP). The design of the NASP ve- 
hicle is a highly complex process. One of the critical tasks involved 
in such a process is the ability of the turbulence model to predict 
hypersonic viscous/inviscid interactions, mixing problems, transition, 
chemical nonequilibria, and a range of other turbulence phenomena. 

Turbulence models are developed on the basis of insight gained from 
experimental and theoretical research. The complexity of turbulence 
requires that the mathematical models be guided by the flow physics in 
a rational and practical approach. Test and validation of new models 
with data of recognized quality is an essential step toward model ac- 
ceptability. A database of high speed turbulent flows is currently being 
completed. This database provides a benchmark for testing and vali- 
dation of new compressible models. Specific topics for the database are 
high-speed attached boundary layers with pressure gradients, super- 
sonic shear layer mixing, and shock wave/boundary layer interactions. 
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The turbulence models will be implemented into the Ames 3-D 
Navier-Stokes codes selected on the basis of several criteria includ- 
ing speed, accuracy, user friendliness, and generality (grid, geome- 
try, boundary conditions). The two-equation models will be validated 
first against the hypersonic database collection. A major consideration 
throughout the research effort is the development of improved com- 
pressibility corrections to the turbulence models and the identification 
of models which are superior in their predictive capabilities. 

Work Statement 

The following is the work statement proposed to be accomplished 
within the first research year: 

1. Evaluate available NASA Ames 3-D RANS codes and make a se- 
lection based on speed, accuracy, user friendliness and generality. 

2. Incorporate candidate compressible turbulent flow models into the 
3-D code and validate against simple attached flows. 

3. Initiate 3-D computations of flows contained in hypersonic database 
collection. 


2 


Work Accomplished 


A summary of the work already accomplished in the first research 
year is outlined below: 

1. Four NASA Ames 3-D RANS codes have been tested and evaluated 
against a flat plate boundary layer flow and an external supersonic 
flow. The codes used in this study are RANS, CNS, FL02, and 
IZTUFF. 

2. Based on speed, accuracy, and user friendliness for this research, 
the code RANS has been selected as the basic code to test the tur- 
bulence models. Work is in progress assisting in the implementation 
of the turbulence models into CNS, FL02 and IZTUFF. 

3. The code RANS has been extended from thin boundary layer to 
full Navier-Stokes. 

4. The k - ui two-equation turbulence model has been implemented 
into the base code. Tests base on supersonic flow at = 5 on a 
cooled flat plate show good agreement with theory and with 2-D 
numerical simulations (Thomas Coakley’s code). 

5. A 24° laminar compression corner flow has been simulated and the 
extend of the separation zone compared against other numerical 
simulations (George Huang, Tom Coakley, Langley’s codes). 

6. A NAS account has been obtained to simulate 3-D turbulent flows. 

7. Work is in progress in the writing of the numerical method of the 
base code including the turbulence model. 
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Code Evaluation Criteria 


The main purpose of the first task was to evaluate and select the 3-D 
NASA-Ames Reynolds Averaged Navier-Stokes (RANS) code or codes 
in order to implement and test compressible turbulence models in the 
simulation of complex high-speed flows with turbulence interactions. 
A guideline of various selection criteria and a guideline of the differ- 
ent numerical methods were designed. The various selection criteria 
included speed, accuracy, user friendliness and generality. The numer- 
ical methods are described including equations, capability, discretiz- 
ing method, convective terms treatment, boundary conditions, viscous 
terms treatment, conservative properties, iteration method, and time 
advance method. 

These guidelines are shown next: 
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Selection Criteria 


We are in the process of selecting the 3-D NASA-Ames Reynolds 
Averaged Navier-Stokes (RANS) code or codes in order to implement 
and test compressible turbulence models in the simulation of complex 
high-speed flows with turbulence interactions. The selection is based 
on various criteria including accuracy, efficiency, user friendliness, and 
generality. 

• Accuracy: 

conservation of mass, momentum, and energy. 

shock capturing capability. 

present validations with experimental data. 

• Efficiency: 

total convergence time to steady-state solutions, 
relative speed: cpu seconds/iteration grid-point, 
robustness, 
stability. 

implicit /explicit numerical scheme, 
implicit/explicit boundary conditions, 
data storage requirements. 
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• User Friendliness: 

user interface. 

instruction manual. 

input parameters definition. 

time-step definition. 

data transfer and management. 

• Generality: 

Navier-Stokes method, PNS method, 
single/multiple grids, 
single/multiple patches in each grid, 
simple/complex geometries, 
simple/complex flows, 
single/general boundary conditions, 
pseudo-time relaxation, time accurate, 
turbulence models. 

chemistry capability (equilibrium, non-equilibrium) . 
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Numerical Methods 


A general description of the numerical method is required including 
the discretizing method, the treatment of the convective terms, bound- 
ary conditions, conservative properties, iteration or relaxation method, 
and time advance method. 

• Equations: 

Mass, momentum, and energy equations. 

Identify Navier-Stokes terms not included. 

Present program status. 

Navier-Stokes, thin layer, Euler, 
elliptic, hyperbolic, parabolic, 
time relaxation, time accurate, 
compressible, incompressible, 
perfect gas, real gas. 
source terms, 
turbulence models. 

• Capability: 

type of problems successfully solved, 
program verifications, 
documentation, user friendliness, 
program availability and support. 
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• Discretizing method: 

finite difference, finite volume, finite element, 
regular or staggered grid, 
simple or multiple grids. 

simple patch or multiple patches in each grid, 
structured, unstructured grids. 

• Convective terms treatment: 

explicit or implicit method. 

central differencing, upwind differencing, or hybrid, 
artificial dissipation, natural dissipation, 
shock treatment, sonic treatment, 
flux-difference splitting, flux-vector splitting, 
first order, higher order. 

TVD, ENO, unconditionally stable. 

• Boundary conditions: 

flow through (inflow, outflow, mixed, fixed), 
wall (adiabatic, fixed temperature, inviscid). 
symmetry (axis, planes, grid-axis), 
characteristic, primitive, or conservative variables, 
explicit or implicit scheme, 
general optional lists or code dependent. 


8 



• Viscous terms treatment: 

no viscous, laminar, and/or turbulent capability, 
explicit or implicit method. 

central differencing, upwind differencing, or hybrid. 

first order, higher order. 

turbulence models incorporated in the code. 

uncoupled or coupled relaxation of turbulence model equations. 

boundary conditions of turbulent quantities. 

• Conservative properties: 

mass, momentum, energy. 

strongly conservative, weakly conservative, non-conservative or 
mass only conservative. 

• Iteration method: 

approximate factorization, relaxation, or hybrid, 
matrix inversion, ADI, DDADI (diagonal-dominant ADI), 
non, over, or under relaxation, 
point, line , plane substitution. 

Newton-Raphson acceleration. 

• Time advance method: 

explicit / implicit . 

first order, higher order. 

pseudo-time relaxation, time accurate. 
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3-D Codes 


The four NASA Ames Research Center codes included in this study 
are RANS, CNS, FL02, and IZTUFF. All these codes have been previ- 
ously used to simulate high speed compressible flows using the three- 
dimensional thin layer Reynolds Averaged Navier-Stokes equations. 
The schemes used in these codes are total variation diminishing (TVD) 
based on flux limiters in the higher-order fluxes. They are spatially 
second or higher-order accurate with first-order accuracy in the vicin- 
ity of shocks. None of the available codes have a manual or a standard 
code description available. These are all research codes under develop- 
ment and they have been previously used to simulate complex flows in 
complex geometries. A brief description of these codes is given below. 

• RANS Code. This code is under development in the Experimen- 
tal Fluid Dynamics Branch by the investigator of this research, 
Dr. Jorge Bardina. It is a 3-D implicit finite difference Reynolds- 
Averaged Navier-Stokes algorithm written in generalized curvilin- 
ear coordinates. It is based on a flux difference splitting method 
which is described in Appendix A under the denomination RANS 
method. It was developed based on the CSCM code 1-4 written by 
this researcher in the Aerothermodynamics Branch coupled with 
ideas obtained from Roe’s flux difference method 5 . During this 
research work, this code has been extended to simulate the full 
Navier-Stokes equations with a general two-equation eddy viscosity 
turbulence model. 

Equations: It solves the 3-D compressible mass, momentum, 
and energy Navier-Stokes equations. Presently, it has the ca- 
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pability to solve the full Navier-Stokes, thin layer, and/or Euler 
equations. 

Capability: It is a research code under development. Its base 
code has been extensively applied in the simulation complex 
3-D hypersonic, supersonic, transonic, and low speed flows on 
complex configurations, internal and external flows. There is 
little code documentation available at present time. The inves- 
tigator in this research has expertise with this code. 
Discretizing method: It uses a finite difference on regular simple 
or multiple-patch grids. 

Convective terms treatment: The convective terms are treated 
using higher order TVD implicit upwind flux difference splitting 
with flux limiters and no explicit artificial dissipation. The 
flux limiters are based on the “minmod” method of Osher and 
Chakravarthy 6 . The averaging procedure is similar to the Roe’s 
scheme, but simpler; it is based on simple arithmetic averaging 
of primitive variables. However, the coefficients of the Jacobian 
matrices are more complex (see Appendix A). It captures shock 
within two to three grid points. 

Boundary conditions: The boundary conditions are implicit and 
they are input through a general optional list of input coef- 
ficients. They include various flow conditions, such as, flow 
through (inflow, outflow, mixed, fixed), wall boundaries (invis- 
cid, viscous, adiabatic, fixed temperature distribution), symme- 
try boundaries (axis, planes), grid axis singularity (in symmetry 
plane) along each coordinate direction. 

Viscous terms treatment: The viscous terms are differenced us- 
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ing an implicit central difference algorithm. They include lam- 
inar and turbulent capability. The turbulence model incorpo- 
rated in this code is the Baldwin and Lomax mixing length 
model 7 . Presently, a two equation k-u turbulence model 8 has 
also been incorporated and tested. 

Conservative properties: Conservation of mass, momentum and 
energy is obtained once convergence is achieved using second or 
higher order differences. 

Iteration method: The iteration procedure is based on an al- 
ternating symmetric plane Gauss-Seidel relaxation procedure 
coupled with a DDADI diagonal-dominant approximate factor- 
ization within the relaxation plane. It provides the capability 
to perform full Navier-Stokes simulations and/or marching pro- 
cedure similar to PNS methods. It includes a Newton-Raphson 
acceleration procedure for faster convergence. It has the capa- 
bility to advance the solution in forward sweeps and/or back- 
ward sweeps. The solution procedure is generally from the in- 
flow plane through the outflow plane firstly, and from the out- 
flow plane through the inflow one secondly. In supersonic zones, 
the best procedure is to advance the solution in the streamwise 
direction only. 

Time advance method: The time advancing is first order with 
capability for uniform time step or local time step based on 
local CFL number. The time step can be increased within the 
boundary layers following T.J. C oakley’s approach to increase 
convergence speed 9 . It can be easily modified to perform time- 
accurate simulations (It has been done and tested previously). 
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• CNS Code. The Compressible Navier-Stokes (CNS) code is under 
development in the Applied Computational Fluids Branch. The 
researchers that provided this code for this research work included 
Dr. Terry Holts, Dr. Thomas Edwards, and Dr. Bradford Bennett. 
It is a 3-D thin layer Navier-Stokes solver in generalized curvilinear 
coordinates on regular zonal grids 10-13 . The flow solver is based on 
F3D flux vector splitting of Steger and Warming 14 , and the zonal 
interface logic is derived from the TNS (transonic Navier-Stokes) 
code 15 . The general treatment of the flux vector splitting is briefly 
described in Appendix A as the Steger and Warming method. 
Equations: It solves the 3-D compressible mass, momentum, 
and energy Euler or thin layer Navier-Stokes equations. It in- 
cludes an air model for equilibrium chemistry. 

Capability: It is a research code under development. It is based 
on the F3D solver that have been extensively applied into the 
simulation of 3-D transonic flows and several hypersonic exter- 
nal flows. There is no manual and little code documentation is 
available. 

Discretizing method: It uses a finite difference on regular zonal 
grids. 

Convective terms treatment: The flux vector splitting of Ste- 
ger and Warming is applied in the streamwise direction, and 
central differencing is used in the crossflow directions. Explicit 
fourth-order (or second-order times pressure gradient ) artificial 
dissipation is added to the right-hand-side in the crossflow di- 
rections to achieve stability with the central difference method. 
Boundary conditions: The boundary conditions are explicit and 
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have to be defined in a special subroutine. Present boundary 
conditions include outflow, freestream, no slip, reflection, axis 
of symmetry, and singularity line at given grid indices. 

Viscous terms treatment: The viscous terms are differenced us- 
ing an implicit central difference algorithm. They include lam- 
inar and turbulent capability. The turbulence model incorpo- 
rated in this code is the Baldwin and Lomax mixing length 
model 7 . 

Conservative properties: Conservation of mass, momentum and 
energy is obtained once convergence is achieved using second or 
higher order differences. 

Iteration method: The flow solver is the F3D which is a two- 
factored flux split algorithm proposed by Ying and Steger 16 . 
The iteration procedure is a LU-ADI method based on the ADI 
approximate factorization of Beam and Warming 17 . The first 
left-hand side operator includes a block lower triangular matrix 
coupled with a block tridiagonal matrix in one crossed flow 
direction. The second left-hand side operator includes a block 
upper triangular matrix coupled with a block tridiagonal matrix 
in the other crossed flow direction. 

Time advance method: The time advancing is first order with 
capability for uniform time step or local time step based on 
zonal CFL number. 
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• FL02 Code. This code is under development in the Aerothermody- 
namics Branch by Dr. Ethiraj Venkatapathy. It is a 3-D implicit, 
finite volume, thin layer, Reynolds Averaged Navier-Stokes method 
in generalized curvilinear coordinates 18-21 . This code is based on 
a LDU-ADI diagonal algorithm developed by Dr. S. Obayashi 18 
with Roe’s averaging 5 and finite differences with MUSCL extrapo- 
lation. Roe’s flux difference splitting method used in this code is 
also briefly described in Appendix A. 

Equations: It solves the 3-D compressible mass, momentum, 
and energy Euler or thin layer Navier-Stokes equations. It 
solves ideal gas and/or chemically frozen flows. 

Capability: It is a research code under development and it has 
been extensively applied into the simulation of 3-D complex 
pl um e flows. There is no manual and little code documentation 
is available. 

Discretizing method: It uses finite differences on regular grids. 
Convective terms treatment: It is a TVD algorithm based on 
Roe’s averaging and upwind flux-difference splitting approxi- 
mate Riemann solver with entropy correction. The MUSCL dif- 
ferencing approach is used to obtain spatially second- or third- 
order accuracy with flux limiters. The flux limiters include the 
“minmod” limiter of Osher and Chakravarthy 6 , second-order 
UNO, second-order differentiable limiter. 

Boundary conditions: The boundary conditions are treated ex- 
plicitly. Present boundary conditions include inflow, outflow, 
freestream, no slip, reflection and axis of symmetry. 

Viscous terms treatment: The viscous terms are differenced us- 
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ing an implicit central difference algorithm. They include lam- 
inar and turbulent capability. The turbulence model incorpo- 
rated in this code is the Baldwin and Lomax mixing length 
model 7 . 

Conservative properties: Conservation of mass, momentum and 
energy is obtained once convergence is achieved using second or 
higher order differences. 

Iteration method: The iteration procedure is based on a scalar- 
matrix operator to obtain steady-state solutions. It is derived 
from an approximate diagonalization of the LU-ADI implicit 
operators in the implicit nonconservative form like Pulliam and 
Chaussee’s diagonal algorithm 22 . It consists of three sweeps 
with a lower, a diagonal, and an upper operator in each coor- 
dinate direction. 

Time advance method: The time advancing is first order with 
capability for uniform time step or local time step. 
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• IZTUFF Code. This code is under develop ment by Dr. Gregory 
Molvik in the Applied Computational Fluids Branch. It is a 3-D 
implicit finite volume thin-layer Reynolds-Averaged Navier-Stokes 
method in generalized curvilinear coordinates 23-24 . This code has 
mainly been used to compute the external forebody flow field of 
3-D bodies, in conjunction with a parabolized Navier-Stokes (PNS) 
code to solve 3-D ideal, equilibrium-, or nonequilibrium-chemistry 
external flows. This code also uses Roe’s method of flux difference 
splitting described in Appendix A. 

Equations: It solves the 3-D compressible mass, momentum, 
and energy thin layer Navier-Stokes equations. 

Capability: It is a research code under development. It has been 
extensively applied in the simulation of complex external con- 
figurations including nonequilibrium chemistry. There is little 
code documentation available at present times. 

Discretizing method: It uses finite volume on regular grids. 
Convective terms treatment: It is a TVD scheme using Roe’s 
upwind-biased flux-difference splitting approximate Riemann 
solver 5 with entropy correction. The “minmod” flux limiter 
of Osher-Chakravarthy 6 is used to obtain spatially second- or 
third-order accuracy. 

Boundary conditions: The boundary conditions are explicit and 
include viscous wall, symmetry plane, and frozen supersonic 
inflow. 

Viscous terms treatment: The viscous terms axe differenced us- 
ing an implicit central difference algorithm. They include lam- 
inar and turbulent capability. The turbulence model incorpo- 
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rated in this code is the Baldwin and Lomax mixing length 
model 7 . 

Conservative properties: Conservation of mass, momentum and 
energy is maintained through the finite volume approach when 
convergence is achieved using second or higher order differences. 
Iteration method: The iteration procedure is based on the 3-D 
ADI approximate factorization of Beam and Warming 17 . 

Time advance method: The time advancing is first order with 
capability for uniform time step or local time step based on 
local CFL number. It can be easily modified to perform time- 
accurate simulations. 
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Codes Performance 


The predictive capabilities of the 3-D Reynolds averaged Navier- 
Stokes codes are well known through the different applications shown 
in available publications (see previous section). The main objective of 
the present work is to analyze and test this codes in order to gain a 
better understanding on their accuracy, efficiency, user friendliness, and 
generality. 

A few tests of inviscid and laminar flows have been simulated. 
These simulations include a supersonic free flow, an hypersonic blunt 
body, and a supersonic flat plate. Some of the results of these simu- 
lations are shown below. In each case, the simulations used the same 
grid and started with the same initial flow field; the flow variables were 
read and rewritten following the particular non-dimensional system of 
each code. 

Figure 1 shows the residual history of a supersonic free flow with 
Moo = 4. In this case, the mesh had llxllxll grid points equally spaced. 
The residuals were defined as the sum of the square residuals of the en- 
ergy equation divided by the number of grid points and the CFL num- 
ber. The energy residual is usually the largest one when it is compared 
with the residuals of the other conservative variables. The CFL factor 
is removed from the residuals in order to get independence from varia- 
tions in the time step or CFL number. IZTUFF achieved a convergence 
of the order of 10 -14 and in about 30 iterations with a CFL = 50; minor 
code modifications were needed in order to achieve convergence. RANS 
achieved a convergence of the order of 10 -n in about 30 iterations with 
a CFL = 50; faster convergence is obtained when the Newton-Raphson 
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acceleration cycle is increased. CNS achieved a convergence of the or- 
der of 10 -11 in about 290 iterations with a TSTEP = 5; simulations with 
higher values of TSTEP produced negative pressure or density values, 
while lower values of this parameter increased the number of iterations 
needed to achieve convergence. FL02 did not achieved a fast conver- 
gence, the residuals were of the order of 7 10 -6 after 9000 iterations with 
a CFL = 50. Some of the differences can be explained by the differences 
in the treatment of the fluxes, however, the results are user dependent 
and code modifications can produce significant improvements. 

Figure 2 shows the residual history of a supersonic laminar flow 
with Moo = 4 and Re = 10 5 . on an adiabatic flat plate. The mesh had 51 
points along the plate, 31 points from the wall through the free stream, 
and 31 points across the main flow. RANS, CNS and FL02 achieved 
a convergence of 10” 9 in about 25, 400, and 3000 iterations, respec- 
tively. RANS was run with CFL = 10000 and four Newton-Raphson 
acceleration cycles. CNS was run with TSTEP — 25 and FL02 was 
run with CFL = 25.0. A slight improvement was observed with FL02 
when CFL = 2, larger or smaller values of CFL decreased the conver- 
gence history of FL02. Dr. E. Venkatapathy achieved convergence in 
about 1000 iterations instead of 3000 iterations when FL02 was run as 
a 2-D plane; this kind of improvement was also observed in the other 
codes. The convergence history of IZTUFF is not included because the 
residuals did not decrease after 1400 iterations with CFL = 0.9; since 
this code is in a developmental stage, it is considered that some code 
modifications are required in this version. 

Although the residual history indicates the number of iterations 
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needed to achieve convergence, different codes and methods require 
different amounts of time and memory to complete each iteration. The 
following table shows the MIPS (millions of instructions per second), 
MFLOPS (millions of floating operations per CPU second), and SPEED 
(CPU micro-seconds per iteration and per grid point) obtained with 
each code in the simulation of the flat plate in the Cray-YMP “eagle” 
at NASA Ames Research Center. 

Table 1. Codes Performance in Cray-YMP 


Code 

MIPS 

MFLOPS 

SPEED 

RANS 

40 

123 

165.6 

CNS 

54 

82 

51.4 

FL02 

34 

124 

50.7 

IZTUFF 

59 

40 

287.1 


This table shows that the CNS and FL02 codes are the fastest codes 
per iteration, the RANS code is about T2 times slower and the IZTUFF 
code is about 5J5 times slower than the other two codes. All the codes 
have been run with a single procesor which has a maximum capability 
of about 300 MFLOPS. This Cray-YMP has 8 processors which can be 
used in parallel. Therefore, the speed of all these codes can be improved 
significantly. 

Previous experience (see previous section) shows that all these codes 
have good shock capturing capabilities and conservation properties. A 
brief comparison of some of the results obtained in the flat plate sim- 
ulations are shown below. Figure 3 shows one of the grid planes along 
the flow direction and the wall boundary. Figures 4a, 4b, 4c and 5a ,5b, 
5c show the pressure contours and the velocity vectors in the vicinity of 
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the inflow obtained with the different codes, respectively. The contours 
and velocity vectors of RANS and FL02 are quite similar, while the 
CNS pressure contours show some oscillations near the wall boundary 
and steeper velocity profiles. Some of these differences are due to the 
different higher-order methods implemented in the codes. 

With regard to user friendliness and generality, all these codes are 
research codes and they are being currently used in the simulation of 
different flows of current scientific and engineering interest. The input 
parameters definitions are explained in each code, although there is 
no instructional manual or user friendly interface. The data transfer 
and management is relatively easy and no major problems were found 
in order to r un these codes when no modifications were required. All 
these codes can simulate complex flows in complex geometries. There 
are strong research efforts undergoing with CNS in the Applied Com- 
putational Fluids Branch and with FL02 in the Aerothermo dynamics 
Branch. The IZTUFF code also includes chemistry capability and a 
two-equation model, however, there is less accumulated experience with 
this code. 
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Turbulence Simulation 


Following the work of Coakley and Huang 8 in the development and 
testing of two-equation turbulence models in 2-D flows, the k-u model 
has been implemented into the 3-D RANS code. Work is in progress 
into the implementation of the k-e baseline models. The model equa- 
tions are well known and they are reported in several publications 8 . 

The simulation of a turbulent adiabatic flat plate flow was chosen 
for evaluation against 2-D results and empirical correlations. The free 
stream Mach number is M 0 0 = 5, the Reynolds numbers are Rei = 
1.4- 10 7 and Re 6 = 10 4 , and the wall temperature is T w /T aw = 1.0. This 
is one of the cases previously selected to analyze turbulence models 8 . 
The numerical simulation was done using 3 parallel planes along the 
streamwise direction. The calculations were done in the middle plane, 
while the lateral planes were used to obtain the Jacobians needed in 
the 3-D code. 

The numerical simulation was run 500 iterations with a CFL = 
500 and required 67 CPU seconds in the Cray-YMP. The speed of the 
calculation was less than 49 CPU micro-seconds per iteration and per 
grid point. This speed is much faster than the one reported in Table 1 
because this was a 2-D calculation, no fluxes were computed and solved 
in the crossed flow direction between the parallel planes. 

Figure 6a shows the 61x61 regular grid used in each parallel plane. 
The grid points are clustered in the wall region and a blowup of the wall 
grid is shown in Figure 6b. Figure 7 shows the normalized temperature 
contours inside the boundary layer. The temperature rises inside the 
boundary layer and is cooled by the cold wall. Figure 8 shows a few 
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selected profiles of velocity vectors. It shows the turbulent boundary 
layer profiles as well as the turning of the velocity vectors due to the 
lip shock. Figure 9 shows the first wall spacing in wall coordinates 
Y + along the wall, it varies between 0.2 and 0.35. Figure 10 shows 
the skin friction C/ distribution along the wall. Figures 11a and lib 
show the turbulent kinetic energy K and mass weighted pK profiles 
against wall distance F+, respectively; and Figures 12a and 12b show 
the specific dissipation rate u> and mass weighted pu profiles against 
wall distance F+, respectively. Figures 13a and 13b show the calcu- 
lated velocity profile compared with empirical correlations, the 1/7 th 
power law and the compressible law of the wall, respectively. Similar 
results are reported in the 2-D calculations 8 . Figure 14 shows the tur- 
bulent eddy viscosity profile with a peak value at about 0.4 boundary 
layer distance, which is inside the range of values between 0.3 and 0.5 
found in other cold- wall boundary layer 8 . Figure 15 shows a relatively 
good agreement between the computed skin friction and the van Driest 
II correlation. Although the agreement is within the observed experi- 
mental uncertainties in cooled-wall flows, the agreement is better in the 
2-D calculations 8 . Finally, Figure 16 shows a relatively good agreement 
between the calculated 3-D and 2-D velocity profiles. This agreement 
is quite good for turbulent calculations done with different codes and 
different numerical methods. 

The results of these simulations support the validation of the tur- 
bulence model incorporated into the 3-D RANS code. 
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Future Work 


The main tasks to be accomplished next are: 

1. Incorporate k — e models into base code, and evaluate performance 
against 2-D flat plate flows and axisymmetric ogive-cylinder-flare , 
Moo = 7, of Kussoy, M.I., et al. 

2. Test and evaluate improved numerics that take into account second 
law of thermodynamics and higher-order TVD flux limiters. 

3. Study and development of compressibility corrections for 3-D tur- 
bulence models: 

3a. Simulate 3-D fin, = 8.2, of Kussoy, M.I. et al. 

3b. Simulate 3-D intersecting SWBLI (2 fins), = 8.3, of Kussoy, 
M.I. et al. 

3c. Simulate 3-D fin, Moo = 6, of Law, C.H. 

3d. Simulate 3-D swept compression corner, Moo = 3, of McKenzie, 
T.M., et al. 

4. Supervise and deliver information needed to implement compress- 
ible two-equation models into the 3-D codes CNS with Brad Ben- 
nett and Tom Edwards, FL02 with Ethiraj Venkatapathy, and 
IZTUFF with Greg Molvik. 
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I. CONSERVATION LAWS 
la. CARTESIAN COORDINATES 

The three-dimensional conservation laws of mass, momentum, and en- 
ergy for the inviseid equations of gasdynamies can he written in the strong 
conservation form as 

dU dFAU ) OFyiU) dFAU) _ (11) 

dt dr dy dz 

where the conservative variables U are 



and the flux vectors F x ,F y ,F z are defined as 

/ pu \ { pv \ / pw \ 

pu 2 + P put ) puw 

F x = puv F,j = pv 2 +p F- = puu» (1.3) 

puw pvu) pw 2 + p 

V u(e + p) / \e(e+p)/ \a’(e+p)/ 

The primitive variables are the density p, the three velocity components 

(u,v,w) and the static pressure p. For a perfect gas, the total energy per 

unit volume, e, is related to the primitive variables according to the equation 

of state which for a perfect gas leads to 

e — p /( 7 - 1) + p(u 2 + v 2 + w 2 )j 2 (1.4) 

where 7 is the ratio of specific heat. 

Ib. GENERAL CURVILINEAR COORDINATES 

The transformation of the conservation laws (eq. 1.1) from (x ,y,z) Carte- 
sian coordinates to general curvilinear coordinates is written in strong 

conservation form as 

dJU dF\{U) dFWJ) dmJ) _ Q 
dt di] dil’ 
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with the flux vectors F, , F 2 . Fi defined as 



and the transformation matrix has the following property 



where I is the identity matrix and ./ is the Jacobian of the coordinate trans- 
formation matrix. 


■J = j'z(y n ~y — yip- if) + //i-- ^ + x <p(yz z it yn~z) 


where the subindex notation of the matrix coefficients means partial differ- 
entiation with respect to the subindex variable. 


Ic. GALILEAN TRANSFORMATION PROPERTY 

Since the Euler equations are invariant under a Galilean transformation, 
it is sufficient to consider 

d JL + Em = o (i.o) 

dt dx 

as the one-dimensional problem associated with the flux-vector splitting and 
the flux-difference vector splitting methods. The one-dimensional conserva- 
tive variables U and flux vector F are defined as 


U = 




F = 


/ pu \ 
I pu 2 +P I 
\u(e + p) J 


( 1 . 10 ) 


In the 1-D case, the primitive variables are the density p, the velocity u and 
the static pressure p. For a perfect gas, the total energy per unit volume, e, 
is related to the primitive variables according to the equation of state as 


e = p/( 7 - 1) + P(“ 2 )/ 2 


( 1 . 11 ) 
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II. FT. r\- VECTOR SPLITTING 


I la. STEGER AND WARMING METHOD 

Steger and Warming 1 have proposed the following flux- vector splitting 
method. Since the flux vector F is an homogeneous function of degree one in 
U. it can be expressed in terms of its Jacobian matrix A as 


where F can be rewritten as 

/ pu \ / pu 

F= pu 2 +p = (3-~i)(pu) 2 /2p + e/h - 1) 

\u(e + p) J \7 e(pu)/p - (7 - l){pii) 3 /2p 2 ) 

or as a function of the conservative variables as 

/ U 2 

F = I (3 — 7 ) C r 2 2 / 2 L r i + Uy/iy — 1) 

\ 7 U3 U2 jU\ — ( 7 — 1 ) G2 3 / 2 Cj 2 ) 


( 2 . 1 ) 


( 2 . 2 ) 


(2.3: 


and the Jacobian matrix .4 becomes 


.4 = 


dF 

W 


0 


-(3 - 7 > 2 / 


2 /9 


1 

(3 - 7 )« 

l)u 2 


-uH + { 7-l)« 3 /2 H-{ 7 

and the total enthalpy per unit mass, H , is 

H = (e+ p)/p = ie/p - (7 - l)u 2 /2 
H = 7p/G(7 - 1) + ^ 2 /2 = c 2 /( 7 - 1) + u 2 /2 



(2.4) 


(2.5a) 

(2.56) 


The flux vector F can be rewritten using a similarity transformation as 

F = AU = SAS~ l U (2.6) 

where A is a diagonal matrix whose coefficients are the eigenvalues of the 
Jacobian matrix A, 


diag A = (u, u + c, u — c) 

c - \frpfp 


12.1 


(2.8: 
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and S is the column-eigenvector matrix 

/ 1 0 0\ /I 1 1 

S = « 10 0 c -c 

\u- 72 U 1/ V° (,2 /(0-l) f'Vb-l) 


/I 

0 

-(7 - !)/T 2 \ 

( i o o\ 


5 _1 = I 0 

1/2 c 

( 7 -l)/2c 2 

-u 1 0 

(2.0b) 

\o 

1 

h— 1 

to 

(7-D/2 c 2 J 

\u 2 / 2 —u 1 / 



The flux vector splitting is defined by the decomposition of the eigenvalue 
matrix A into a non-negative and a non-positive diagonal matrix. 



A = A+ + 

A“ 

(2.10a 

+ = (A+ 

1 A |)/2 

(2.10/) 

" = (A- 

1 A |)/2 

(2.10c 


where | A | is defined as a matrix whose coefficients are equal to the absolute 
value of the corresponding coefficients of the matrix A. Thus, A + includes 
only the non-negative eigenvalues of A and A~ includes only the non-positive 
eigenvalues of A, respectively. Following this splitting method, the flux vector 
F is rewritten as 


F = S( A+ + A ~)S~ X U = (.4+ + A~)U = F + + F~ (2.11) 


The spatial derivatives of F + and F~ are usually approximated with the 
standard implicit first- and second-order backward and forward difference 
operators, respectively. Flux limiters are introduced in the higher-order dif- 
ference operators to eliminate oscillations in the presence of discontinuities, 
such as shock waves. The Euler equations (1.8) become 


dU dF + dF~ 

dt + dx ^ dx 


( 2 . 12 ) 


and the difference operators are 


KK = (Ft - F+.ySx + 4,- (Ft - 2F+ ! + Ft_ t )/26x 
6+F- = (F~ +1 - F-)/6x - 4>f(F- - 2 F~ +l + F~ +2 )/26x 
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(2.13) 

(2.14) 



with first-order when o = 0 and second-order when o = 1. Other higher- 
order difference operators are expressed in a similar form, such as Fromm 
scheme and third-order biased scheme. Equations (2.13) and (2.14) can also 
Ire written as a extrapolation at cell interfaces 

^F- = (F- 1 / 2 -F- 1/2 )/^.r 


(2.15) 

(2.1G) 


where 




(2.17) 

(2.18) 


with the appropriate definition of <j> at cell interfaces. 

An alternative approach proposed by van Leer 2 is the MI SC L method. 
In this method, the variables U are extrapolated toward the interfaces and 
then the flux vectors F + and F~ are evaluated. Therefore 


6 Z F, = ( F+(F,: 


+i/ 2 ' 


F + (U-_ l/2 ) + F 


-(Ut 


t+l / 2 


)~F-{UU 


/ 2 > 


)/S.r 


(2.19) 


and 


U ti/2 = Vi + - Ui )/ 2 (2.20) 

Ur +1/2 = U i + 4>-(U i -U i . 1 )/2 (2.21) 
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III,. VAX LEER METHOD 


Although the flux vector F is a continuous differentiable function, the 
fluxes split F + and F~ of Steger and Wanning described in the above previous 
section are not differentiable at zeros of the eigenvalues. A method proposed 
by van Leer defines the split fluxes as continuous differentiable functions 
including sonic and stagnation points. The flux vector F expressed as a 
function of density p. sound speed r, and Mach number M = u/c becomes 

/ pcM 

F= pc 1 ( M 2 + 1 / 7 ) 

\pc 3 M ( M 2 / 2 + 1 /( 7 - 1 ) 

The formulas for the flux split given in terms of the Mach number are for 
supersonic flows, | M |> 1, 

F + = F F~ = 0 for M > 1 

F + = 0 F~ = F for M < -1 

and for subsonic flow. | M |< 1, the flux vectors are defined a 

( ±pc((M ± l)/2) 2 

F ± = ±pc 2 ((M ± l)/2) 2 ((7 - 1 )M ± 2)h 

\±pc 3 ((M ± 1)/2) 2 ((7 - 1)M ± 2) 2 /2( 7 2 - 1) 


(2.23) 

(2.24) 
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III. FLUX-DIFFERENCE SPLITTING 


The one-dimensional conservation laws of mass, momentum, and energy 
for the inviscid gasdynamics equations can be written in a discrete space as 

6,U = -At 6 t F{U)/Ax (3.1) 


There are significant differences in the numerical methods between most 
finite volume and finite difference approaches. One of these differences is the 
location of the cell areas. 

• Finite Volume method. 

In a finite volume approach, the flux difference splitting are defined as 
follows 

S x Fi = 6 X F+ + S X F~ (3.2) 

where F, = F(U,) and the flux differences are 


UU + =F+ l -F+ l 


l+i 


and 


6,F-=F- -F~_ 


(3.3) 


and 


Ax = (x.+i - Xj_ i )/2 


(3-4) 


The flux difference may be rewritten also as 


SF t = h i+ i 


F 


(3.5) 


• Finite Difference method. 

In a finite difference approach, the flux difference splitting are defined as 


S X F = 6 X F+ | S X F~ 
A.r A.r + Ax - 


(3.6) 
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and a first-order upwind method gives the following equations 


( 7 F+ = f+- f+, 

(3.7 

~ j — F i 

(3.S 

= .r, - x,_i 

(3.9 

A r = .r,+ 1 — x t 

(3.10 


In both methods, the fluxes are split in delta form based on a Jacobian 
matrix A which follows ‘Roe's property U.' 

6F = AbU (3.11) 

Under a similarity transformation, the matrix A is expressed as 

A = S AS' 1 (3.12) 

where A is a diagonal matrix with the eigenvalues of the matrix A as diagonal 
coefficients, the matrix S is composed by the eigenvectors of the matrix A. 
The matrix S~ l 8U is also the expression of the characteristic variable rep- 
resentation of the conservative variables U. The positive and negative flux 
differences are defined according to the sign of the eigenvalues as 

6F? = A + (Ui ~ Ui- 1 ) and 6F~ = A~ (U t+i - U,) (3.13) 

where 

A + = S A + S~ l and A~ = S A - S 1 (3.14) 

A + = (A+ | A |)/2 and A - = (A- | A |)/2 (3.15) 

In both methods, finite volume and finite difference, the flux differences 
may be combined into a equivalent expression as (H.C. Yee, NASA TM 
101088) 

SF l — hi , i — hi i 
i-r 2 i 2 


(3.19) 


whore 


(3.20) 


= ( -F’i + 1 + Fi — | .4 | (l t + 1 — 1 1 ))/'- 

where | .4 | is the Jacobian matrix with the absolute value of its eigenvalue 
coefficients, 

| .4 |= 5 | A | S -1 (3.21) 

The main differences between the numerical methods are due to the choice 
of finite volume or finite difference approach, the definition of the Jacobian 
matrix ,4, and the implicit operator procedures.. 
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Ilia. ROE METHOD 


The Roe's method described here is a finite volume, flux difference split- 
ting method using Roe s averaged variables. This method uses the above 
defined equations (3.19, 3.20. and 3.21) and it is based on the property that 
the the flux F is an homogeneous function of degree one of the conservative 
variables U. Therefore 

F = AU (3«.l) 

where A is the Jacobian matrix of the flux vector F. Thus, at each nodal 
point i 

F t = F(Ui) = A,U, (Zo.2) 

On the other hand, at the cell interfaces between the cells r, and ,r l+1 . Roe 
defined a set of averaged variables to define the coefficients of the Jacobian 
matrix .4 subject to the following condition 

SF = Fi+\ -Fi = A(U,+\ - Ui) = A 6U = SAS~'6U (3c. 3) 

The average density is defined as the geometric-averaged value of the nodal 
points, the averaged velocity components and total enthalpy are defined as 
weight ed-averaged nodal values 


* + i - \/Pi+iPi 

( 3o.4) 

* + i = (Wu l+1 + Ui )/(W + 1) 

(3a.o) 

t * + . =(WH l+l +H t )/(W + l) 

(3c. 6) 

+ 

11 

♦ 

1 

C 

* 

ro 

to 

(3c. 7) 

W = yj pi+ij pi 

(3c. 8) 


Thus, the Jacobian matrix becomes 


0 

-(3- 7 )u* 2 /2 
-u*H* + (7 — 1 ) u * 3 /2 


(3 
H* - 


0 

h-i: 
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and the column-eigenvector matrix and its inverse are 



and the diagonal matrix 


o 0W1 i i 

10 0 c -r 

u* 1/ \0 r*' 2 /h-l) ^ 2 /(0-l) 


0 

-(7-D/<\ 

( 1 

0 

1/2 c* 

h-l)/2c' 2 


1 

-1/2 c* 

( 7 -l)/2 c* 2 / 

\ u* 2 /2 

— u* 

A is 






(3a. 10) 


(3a.ll) 


diag A = (a*, u* + c* , u* — c*) 


(3a. 12) 
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III1). CSCM METHOD 


This method is a finite difference, flux difference splitting method us- 
ing 'Roe's property U'. This method uses the above defined equations (3.6 
through 3.18) and it is based on the following property 

6F = A6V (36.1 ) 


The flux difference splitting is defined based on 


6F ± = 5/ ± 5 _ 1 .4 6U ( 36.2 ) 

where I is the identity matrix and it is splitted into a non-negative matrix 
and a non-positive matrix according to the sign of the coefficients 


diag / ± = ±1 • sign(u, u + c, u — c) 


(36.3) 


where the overbar denotes arithmetic averaging. 
The matrix .4 is defined as 


/ 


A 


0 


(7 - 1)(« 2 


\u(7( 


u 2 - 


sl)-u£2L 
2 ' p 

bK — tp 




p( 7 ~* 1 ) 


(2-7)« + f 

H _ ,7 (~ r . — JE ) j H — 

2 U( ^ U p’^pi 7 - 1 ) 


0 

7 - 1 
7“ 


(36.4) 


and the matrix S represents a transformation matrix from the conservative 
variables to characteristic variables, 


5 = 



S~ l = 



00 \ f-p m p/ 2 \ 

1 0 0 pci 2 -pet 2 

u 1/ \ 0 7p/2(7-l) tT/2(7 — 1) / 

0 (7 - 1)/7P\ / 1 0 0\ 

1/pc (7 - 1)/7P J ( 1 0 

-l/pc (7 — 1)/7P / \ u 2 — u 2 /2 -u 1/ 


(36.5) 


( 36 . 6 ) 


61 


III. . RAXS METHOD 

The RAXS's method is similar to Roe's dux difference splitting method, 
however it is based on simple arithmetic averaging of primitive variables. The 
flux differences are also defined as 

6F = F i+ 1 -Fi = A{U i+ 1 - Ui) = A6U = SAS~'6U (3c.l) 

in order to preserve the conservative property of the fluxes. The column- 
eigenvector matrix S and its inverse are 

/ 1 0 0\ (-p p/2 p/2 _ \ 

5 = I _u 1 Olio p(c-d)/ 2 -p(c + d)/ 2 1 (3/'.2) 

\m 2 /2 a l/\0 tP/2( 7 -1) 7p/2(7 - 1 )/ 

/-1/p o (7 — i)/tP \ / 1 0 0\ 

5- J = 0 l/pc (l+d/c)(7-l)/7P 10 (3r ‘ 3) 

\ 0 -1 /pc (1 - d/c )(7 - l)hpj \ U 2 -U 2 / 2 -u 1/ 

and the diagonal matrix A is defined with the eigenvalues of the matrix A as 

diagA = (u, u + c, u — c) (3c.4) 

and the coefficients 

(3e.5) 
(3c. 6) 
(3c. 7) 


ti = (u + pu/p)/ 2 
d = (u - pu/p)/ 2 
c = \mP/p) + d 2 
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RANS METHOD 


I. CONSERVATION LAWS 

Tli*' conservation laws of mass, momentum, and energy for the 3-D 
Navier-Stokes equations can be expressed as 

dU 0F\ OF 2 OF* n M 

1 H -\ =0 la 

9t dr, d.r 2 a.r 3 

or in compressed vector notation (where repeated sub-indices in any term 
imply summation over the index range, sub-index t following a comma imply 
partial differentiation with respect to time, and sub-index numbers following 
a comma imply partial differentiation with respect to the respective spatial 
coordinate directions) as 

U,t + Fj,j= 0 (ll>) 

The vector U represents the conservative dependent variables, 

U = (p,pu i ,pu 2 ,pu 3 ,e) ('-) 

and Fj are the flux vectors in the respective Cartesian coordinates direction 

Fj =(pu 3 , pu\Uj 4- phj — V] j,pu 2 Uj + pS 2j — r- 2 j, 

pu 3 uj + pS 3 j - r :i j , (e + p)uj - u,Tij - kTj) (3) 

p is the fluid density, u } are the velocity components in each coordinate di- 
rection, p is the static pressure, k is the thermal conductivity, T is the fluid 
temperature, and e is the total energy. 

The viscous stress tensor r tJ including the turbulent eddy viscosity is 
defined as 

Tij = 2(p + Pr)(Sij - SijSnn/ 3) (4) 

where S tJ is the second-order isotropic Kronecker delta. 
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For a perfect gas. the pressure' is related to the total energy as 


( = p/il ~ 1 ) 4 - pu.jUj/2 ( 5 ) 

where is the ratio of the specific heats at constant pressure and volume, and 
the temperature is related to the density and pressure through the equation 
of state 

p = RpT (G) 

The conservation laws written in generalized curvilinear coordinates = 
bUj) ;ir(1 

U,i + £1 j Fjj = 0 * ~ * 

where i represents the curvilinear coordinate direction £,. The conservation 
laws can be rewritten in weakly-conservative form as 

(JU),t+7^~jF jrt = 0 

or in strongly-conservation form as 

where J is the Jacobian of the coordinate transformation 

J = *ijk x \~i x 2~i X Z~k 

and e l}k is the third-order isotropic alternating tensor. The metric coefficients 
are evaluated at the interface between the grid points of the spatial difference, 
according to their mathematical definition as 

J £i,s = 0-G f ijk fjnm X n,j X m,k ' " ^ 

The present version of the RANS code uses the weakly-conservative form 
shown in the above equation 8 . Different tests and experiments in subsonic, 


(S) 


( 9 ) 


65 



transonic, and supersonic flows have shown that the conservation properties 
are maintained once convergence is achieved with a second- or higher-order 
scheme. An strongly-conservative version of the code is in progress. 

II. FLUX-DIFFERENCES 

All flux differences are treated implicitly in order to increase stability 
and to be able to use large increments of time or CFL numbers. The numer- 
ical scheme for the viscous fluxes is second-order central difference, while the 
numerical scheme for the inviscid fluxes is higher-order T\ D upwind flux- 
difference splitting. The higher-order T\ D scheme is based on the Oshei- 
Chakravarthy ” minmof method applied to the flux differences; it has the 
capability of to represent various higher-order differences: first-order up- 
wind, second-order upwind, third-order upwind biased, second-order Fromm 
scheme, and other combinations of upwind and central differences. 

The differentiation of the inviscid fluxes is similar to Roe’s flux differ- 
ence splitting method, however it is based on simple arithmetic averaging 
of primitive variables. Each flux difference term iy t of the conservation law 
equations is defined as 

SF = F 1+i - F- = A(T- +1 - U-) = A6U (12) 

in order to preserve the conservative property of the fluxes (the subindex j is 
omitted in equation 12 for simplicity). The Jacobian matrix A is decomposed 
using a similarity transformation as 

A = SAS -1 (13) 

where A is a diagonal matrix whose coefficients are the eigenvalues of the 
matrix A, and 5 is the column-eigenvector matrix of A. The matrix 5 is 
further decomposed in two matrices 

S = MT (14) 
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where M represents the transformation matrix of conservative and primi- 
tive variable differences, and T represents the transformation matrix between 
characteristic and primitive variable differences. For the flux differences in 
the coordinate direction, these matrices are defined as 
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where the arithmetic averaging is denoted by the overbar symbol. Although 
the density p could be factored in the T and T~' transformation matrices, as it 
has been pointed out by T..J. Coakley, the density factor has been introduced 
in order to write the characteristic variables in nondimensional form. The 
prime symbol (x') is used to denote unity magnitude of the vector, thus 

=^,5/ 161 

= /I6I 

16 1 = \J + 66 + (j ,3 
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( 17 a) 

( 176 ) 

( 17 c) 


with no summation of repeated indices in the above equation (17). The 
hat symbol (./•) is used to denote the inverse matrix coefficients of the noil- 
orthogonal matrices. 

■b,2 - ^,3(^.2 •^,3))(I^|I6I/./|6I) 2 (1S«) 

— ( .c 1 2 — 1 j( ,r i,'2 ■ )) ( 1 — (-C, '2 ' - r i,3 1 ) ( ISI>) 

= (^,3?.,3 -^/,,2)(l6H6l/-/|ed) (13C) 

which are all equivalent expressions. 

The diagonal coefficients of the eigenvalue matrix A for the flux differences 
in the coordinate direction are 

diagA = hjuj, J uj uj, 0.5 + Uj) H- c. 0.5£, Jiij + u } ) - r) (19) 

where 

c = J 2 (20) 

d = l hj 6u } ( 21 ) 

Suj = ( uj - puj/p )/ 2 ( 22 ) 

Uj = (Cl] +pu J /p)/2 ( 23 ) 

The inviscid first-order flux differences terms (see equation 12) are split 
according to the sign of the eigenvalues as 


6F = VF + + AF~ (24a) 

VF+ = .4+ VU = S\ + S~' VU (24 b) 

A F~ = A~ A U = 5A _ 5 _1 A U (24c) 

where V is the backward-difference operator and A is the forward-difference 
operator. The coefficients of the diagonal matrices A ± are defined as 

+ |A|)/2 (25a) 

A" = (A-|A|)/2 (256) 
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Unconditionally stable Euler implicit methods are constructed by forward 
upwind differencing flux differences with positive eigenvalues and backward 
differencing flux differences with negative eigenvalues. 

III. HIGHER-ORDER TVD FLUXES 

Higlier-order spatial TVD flux differences in the right-hand-side of the 
inviscid terms of the conservation law equations are defined by using the 
” minmod" limiter of Osher and Chakravarthy applied on the complete split 
flux difference instead of the characteristic variable difference. For a given 
direction and at the j grid point. 

VF = A,_,F+ + A J+ ,F- ^ 

+ " [ + br ( ( A ' + F~ - A m f ~ > + <■ Vi- F+ - A V f+ > ) 

+ by (*‘V7 F+ - VjC) +(A,Wf- - A j+i F-)) } (20) 

where the first-order terms, the second-order upwind-difference increment 
terms, and the second-order central difference increment terms are shown in 
the first., second, and third line of the above equation 26, respectively. 0 < 
c < 1 is an input factor that allows to change from first-order to higher-order 
differences monotonically. A F and A F are the limited split flux differences 
defined by 

A ;+ iF + = rninmod( A J+ ^ F + , /?A j_^F + ) 

A ;+ |F _ = minmod( A J+ | F~ , (3 A }+ iF~) 

A-sF+j = minmod(A-3F + , {3A._\ F + ) 

J 2 J 2 J 2 

= 77iinmod(A i _±F~, (3 A ■ . if") (27) 

J 2 / v J 2 j ' 2 

The minmod function is defined as 

minmod{a,b ) = sign(a) jnax(0, min( |a|, bsign(a))) (28) 
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and iv turns a value of zero if the product ah < 0. or a value equal to the 
argument with the minimum magnitude if both arguments a and 6 have the 
same sign. 3 is a compression parameter determined in the range given by 

1 < 0 < - ° (29) 

- “ 1-0 

The maximum value of 3 is used because it introduces less dissipation. The 
numerical scheme is second-order upwind when 0 = -1 (3 = 2). Fromm's 
scheme when 0 = 0 (3 = 3), third-order upwind biased when 0 - 1/3 (3 = 4). 
and second-order central-difference when 0 = 1 (3 —> oo). 

For example, the second-order flux-difference split upwind scheme with 
c = 1 is 

+lA 1+i F- - \ a^F- (30) 

where the minmod function 

0 < 0.3 F+ = minmod(A. 3 f+ 2A,_.F + ) < 2A. i F + (31a) 

— J 2 J 2 J 2 J 2 

0 > i\j + 3F~=minmod(Aj + sF~,2Aj + iF~) > 2A J+ iF (316) 

Thus, the TVD scheme limits the split fluxes with a factor between 0.5 and 
1.5 times the first-order flux-difference, 
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-A i+ xF~ 
2 7+2 

> -A,i F~ 
— 2 J+ 2 

- -A j+ 3 F- 
2 ^ 2 

> -A, , i F~ 
~ 2 J+2 


IV. IMPLICIT NUMERICAL METHOD 

The numerical method is based on an implicit ” method of planes” sym- 
metric Gauss-Seidell relaxation scheme. The data is conveniently stored on 
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successive planes along the streamwise coordinate’, and the system of equa- 
tions is solved in each successive plane along the forward direction, first, 
and along the backward direction, afterwards. In each plane, the solution is 
obtained by using a two level pseudo time dependent relaxation procedure 
based on a diagonally dominant approximate factorization DDADI. The space 
marching alternating directional sweeps in the streamwise direction are von 
Neumann unconditionally stable for zones of subsonic and streamwise sep- 
arated and reversed flows as well as supersonic flow. The space marching 
method results in improved propagation of nonlinear effects to accelerate 
convergence to steady state, much as do the more restrictive PNS techniques. 

The diagonal dominant approximate factorization of the left-hand-side of 
the conservation law equations including the implicit viscous terms leads to 
the following block tridiagonal equation sequence for the <h plane relaxation 
method 


(-At , D, A7 )6U* = -RHS n ' n+1 (33 a) 

(~A^D,A- a )6U = -D6U* (33 b) 

The diagonally dominant matrix D involves the first-order split Jacobian 
matrices and the Jacobian matrices of the viscous terms of each coordinate 
direction 

D = I + 4 - y, + y - V, + A* b - A- (34) 

and the solution is updated from time step n to time step n + 1 

U n+l =U n +6U n (35) 

Observe that the RHS of equation 33a has an exponent n,n + 1 because some 
terms in the streamwise direction are already updated at time step n + 1 due 
to the plane relaxation procedure. 
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A Xewton-Raphson acceleration procedure is obtained by solving each 
plane twice or more times in each relaxation sequence. This procedure pro- 
duces significant improvement in the propagation of nonlinear effects due to 
the nonlinearity of the .Jacobian matrix .4. 
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